Effects of Pauli paramagnetism on superconducting vortex phase diagram in strong 

fields 



cn 

O 

o 

(N 
^■ 



o 
o 

u 

Oh' 



X3 
O 

> 
O 

O 

o 
m 
o 



I 

o 
o 



X 
J3 



Hiroto Adachi and Ryusuke Ikeda 

Department of Physics, Kyoto University, Kyoto 606-8502, 
(Dated: February 2, 2008) 



Japan 



The Ginzburg-Landau (GL) functional and the resuhant phase diagram of quasi two-dimensional 
superconductors in strong fields perpendicular to the layers, where the Pauli paramagnetic depairing 
is not negligible, are examined in details by assuming the weak coupling BCS model with a d^2_y2- 
like pairing. It is found that the temperature at which the mean field (MF) transition at Hc2 changes 
into a discontinuous one lies much above another temperature at which a line -ffpFLO (T) of transition 
to a helical FFLO-like vortex solid may begin to appear. In addition to MF results, details of a real 
phase diagram near _ffc2(T')-line are examined based on a theoretical argument and Monte Carlo 
simulations, and it is found that the MF discontinuous transition is changed due to the fluctuation 
into a crossover which is nearly discontinuous in systems with weak enough fluctuation. These 
results are consistent both with the MF discontinuous behavior and a suggestion of -ffFFLo(r) in 
the heavy fermion superconductor CeCoIus with weak fluctuation and with their absence in cuprate 
and organic superconductors with strong fluctuation. 



PACS numbers: 74.25.Dw, 74.40.+k, 74.70.Tx 



I. INTRODUCTION 

Traditionally, effects of Pauli paramagnetism on su- 
perconductors with a spin-singlet Cooper-pairing have 
been discussed by simply focusing on two energy scales 
the superconducting (SC) condensation energy and 
the Zeeman energy preventing the singlet pairing. This 
is a reasonable explanation on the first order transition 
(FOT) in an exceptional case with no orbital depairing 
creating field- induced vortices, i.e., a thin- film supercon- 
ductor in parallel fields Q. Further, there is also a possi- 
bility of a structural transition within the Meissner phase 
into the so-called FFLO state 0, 0] with a periodic mod- 
ulation, induced by the spin (paramagnetic) depairing, 
of the SC order parameter. However, the orbital depair- 
ing effects, i.e., the field-induced vortices, are inevitably 
present in most of cases of a bulk type II superconductor 
under a strong field, including a layered material under a 
field parallel to the superconducting layers 7]. Hence, we 
encounter quite a different issue from that in the works 
Q IM Q IE 111; that is, effects of the spin depairing on 
the vortex state which has no Meissner response. Since 
the number of vortices is determined only by the mag- 
netic field strength and system sizes, treating the orbital 
depairing as a perturbation for the case with no orbital 
depairing is not valid. 

At present, it is well understood [E @] that, in lower 
fields where the spin depairing is negligible, the H-T 
phase diagram for the vortex states is drastically changed 
by including the SC fiuctuation neglected in the mean 
field (MF) approximation. A typical one among such 
drastic fluctuation effects is the fact that the second or- 
der MF transition at Hc2 is not realized as a consequence 
of the fluctuation and gives way to a weak first order 
transition lying below Hc2 between the vortex solid and 



the vortex liquid region which needs not be distinguished 
from the normal state. At least theoretically, it is impor- 
tant to extend this issue to the strong field region in 
which the spin depairing is not negligible and the MF 
transition at Hc2 may be discontinuous. 

Through previous MF works on the vortex states of su- 
perconductors with paramagnetic depairing [i(tI llll Il2l | , 
however, one notices that even the H-T phase diagram 
in the MF approximation is an unsettled issue. For in- 
stance, a FOT induced by the spin depairing was ex- 
pected through a calculation in dirty limit 12], to the 
best of our knowledge, contrary to the experimental facts. 
Further, even in pure (clean) limit, there are no consen- 
sus on the MF phase diagram. In strong fields of our 
interest, any Meissner phase (i.e., any phase occurring 
with no orbital depairing) is not possible, and we ex- 
pect just some vortex solids, such as the ordinary solid 
consisting of straight vortex lines and an, if any, FFLO- 
like solid state with a periodic modulation along the ap- 
plied field, as SC ground states in clean limit with no 
defects leading to a vortex-pinning. Let us call a transi- 
tion curve between the above-mentioned two vortex solid 
states as i/FFLo(2^)- In ref.^3|) the transition at TJpFLO 
was argued to be of first order with no detailed calcu- 
lation, while it was obtained as a second order one in 
ref.jll] where the orbital depairing represented by the 
gauge-invariant gradient is treated perturbatively. Fur- 
ther, the temperature T* at which the MF transition 
at Hr7 changes into a discontinuous one was concluded 
there [ijj to lie much below another temperature Tfflo 
at which Hppi^o(T) and Hc2{T) branch. In addition, 
the i/pFLO (r)-line is often suggested to be insensitive to 
T. We note that all of conflicting results raised above 
were obtained in the same model, i.e., the simplest weak- 
coupling BCS model. 



2 



In the present paper, we consider the SC phase di- 
agram of quasi two-dimensional (2D) superconductors 
with a singlet d-wave pairing under a high field perpen- 
dicular to the SC layers where both the orbital- and spin- 
depairings are not negligible. In Sec. II, a MF calculation 
for the weak-coupling BCS model is first performed by 
focusing primarily on the region where the SC order pa- 
rameter is described within the lowest Landau level (LLL 
or = LL). We demonstrate by treating the two de- 
pairing effects on the same footing that most of the pre- 
vious MF results mentioned above may be significantly 
changed. Our analysis is different from that in ref.jnl in 
that the gauge-invariant gradient is fully incorporated. 
This nonperturbative inclusion of the orbital depairing 
should become important upon cooling and with increas- 
ing field, although, instead, multiple numerical integrals 
have to be performed to accomplish it. Further, bearing 
a comparison with data in real systems in mind, a weak 
impurity scattering should be incorporated because the 
HpFhoiT) line which may appear in high H and low T 
is found to be quite sensitive to the impurity strength. 
Consequently, the MF phase diagram is determined by a 
competition among the two depairing effects and the im- 
purity strength. We find that the temperature T* always 
lies above other two characteristic temperatures Tfflo 
and Tnext below which the vortex solid will be described 
by the next lowest Landau level {N = 1 LL). In addition 
to these studies of the MF phase diagram, linear response 
and elastic properties in a possible helical FFLO-like vor- 
tex solid will be examined to give a result useful for search 
of such a helical solid state. 

Next, in Sec. Ill, we address effects of the SC fluctua- 
tion on the phase diagram by neglecting vortex pinning 
effects. We notice that a theoretical consideration on the 
fluctuation effects on the orderings in vortex states below 
Hc2 is unaffected by a change of the GL model brought by 
the presence of the spin depairing and naturally leads to 
the conjecture that the FOT obtained in the MF approxi- 
mation (MF-FOT) at Hc2 will not occur as a true FOT in 
real systems. Actually, otherwise, the high field portion 
of the H-T phase diagram would not become compatible 
with its low field portion in which the absence of the sec- 
ond order MF transition at Hc2 is well-established 0, . 
The only transition in the case with no pinning disorder 
should be the melting transition of a vortex solid. In the 
present case with a MF-FOT, however, the GL model 
needs to have higher order (nonlinear) terms other than 
the quartic term, and this fact makes an analytic study 
more involved. For this reason, we have chosen to per- 
form Monte Carlo simulations on a GL model justified 
through our microscopic analysis in order to examine the 
true phase diagram. Our simulation results are consis- 
tent with the above-mentioned conjecture that the MF- 
FOT at Hc2 is changed due to the fluctuation into a 
crossover of which the width is narrower in systems with 
weaker fluctuation. Together with the impurity-induced 



disappearance of the MF-FOT, this fluctuation- induced 
broadening of the sharp behavior reflecting the MF-FOT 
explains why the nearly discontinuous behavior at Hc2 
has not been observed so far in, except the recent ob- 
servations in CeCoIns [H H [H E [13 , most of bulk 
type II superconductors with a spin-singlet pairing. On 
the other hand, we find that, in spite of the absence of 
the genuine FOT at i?c2, a hysteresis may appear near 
Hc2 in simulations for cases with weaker fluctuation (or 
equivalently, at very low temperatures) as a result of an 
incomplete relaxation at finite time scales. 

Finally in Sec. IV, experimental facts in CeCoIns un- 
der H II c and other materials are discussed based on the 
results obtained in Sec. II and III. It is conjectured there 
that the SC fluctuation makes the HFFi^o{T)-\me and 
the FOT-like behavior at Hc2 unobservable in strongly 
fluctuating systems such as organic materials. Further, 
it is also pointed out there that, although microscopic 
results in Sec. II are not applicable to the H _L c case in 
CeCoIns, recent specific heat data 0,0| can be under- 
stood within our results in Sec. III. 

We use the unit h — c = kB = 1 throughout the 
manuscript. 



II. MEAN FIELD ANALYSIS 

Following Klemm et a/.|0|, we start from a BCS hamil- 
tonian for a quasi 2D superconductor under a nonzero 
magnetic field perpendicular to the basal plane 



Ti. — Tin + J 



2 ^ 



Sjt(q^)BJ(q^).(l) 



Here, g is the attractive interaction strength, ip'j{r±) = 
(l/y/fls) X]p^ aJ(p^)e'P^ ''^ is the annihilation operator 
of a quasi-particle with spin a (= +1 or —1) at the in- 
plane position rj^ on the ^'-th plane, and s and O are 
the interlayer spacing and the area of a layer, respec- 
tively. The pair-field operator is defined by BJ{c[±) — 
WpaJ'^(— p2)aj(pj|), where Wp is the orbital part 
of the pairing-function and, in the case of da;2_j,2-pairing, 
is written as ^/2{p^ — Py), where p is the unit vector in 
the p^ direction, and p^ implies pj^±q/2. The first 
term Tig in eq.(^ represents in-plane motions of quasi- 
particles. 



Ho - sY^J d?r^{^'^\v^) 



(-iV^ + |e|A) 
2m 



■<(rj 



+^f(r^)^«,(r^)-a/J^J(rj.)J, (2) 

where the gauge field A has no out-of-plane component, 
/ = ^.qH is Zeeman energy, and m is the effective mass 
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of a quasi-particle. The strength of the paramagnetic de- 
pairing is measured by fj,oH°2^ (0) /2ttTco corresponding 
to the Maki parameter except a numerical factor, where 
H°2^ (T) is the MF transition curve in the case neglecting 
both the paramagnetic depairing and impurity effects, 
and Tco is the zero field {H = 0) transition temperature 
in the pure limit. Throughout the numerical calculations 
in this and next sections, we will choose iJoH°^'"{0)/2TrTco 
to be 0.8 by bearing a comparison with observations 
in CeCoIns in mind (see sec. IV). The random potential 
Vj{r±) obeys the following Gaussian ensemble: Vj{r±) = 
0, Vj{r±)vj>{r'^) = S{rj_ - r'^)6jj> /2'kN{Q)t, where iV(0) 
is a 2D density of state at the Fermi surface, and the 
elastic scattering rate. The second term represents the 
inter-plane hopping, 

= IdVx^^^f (r±VJ+i(r^)+H.c.y(3) 

where H.c. denotes Hermitian conjugate. We use 
the familiar quasi-classical approximation for the single- 
particle propagator 



G" {v,T ) — G" (r — r )e -I" 



(4) 



Here, the Green's function in zero field, °(r), is given 
as the Fourier transform of the expression 



GeAp) = ie^ - (Cpi + J cos(p^s)) 



(5) 



where e denotes the Matsubara frequency £„ = 2'KT{n + 
1/2) in this paper, = e+Se/ {2T)—iaI , and = sgn(£). 
Since we take account of the paramagnetic depairing sup- 
pressing the MF upper critical field, the use of the quasi- 
classical approximation, valid if kpru 2> 1, needs not to 
be questioned, where fci? is the Fermi momentum, and 
th = (2|e|i7)~^/^ is the magnetic length. According to 
the familiar Stratonovich-Hubbard procedure to intro- 
duce the pair- field A, we can construct a GL functional. 
Throughout this paper, we neglect the internal gauge 
fluctuation and examine the resulting GL model in type 
II limit. The validity of this approximation will be ex- 
plained in discussing the MF phase diagram later. 



A. Quadratic term 

The quadratic term of the GL functional is given by 
•^2 = E/ rf'^xA;(r^)(^i|~i^2(n,g,)) A,^(r^), 

- . . (6) 
where Aj = (iViaycr) "^^^ Ag^e"'^"-' denotes the pair- 
field on the j-th SC plane, Mayer is the number of SC 




0.4 0.6 

T/Tco 

FIG. 1: Comparison between iifc2(r)-curves for = 
0.1 determined from ao(0) = by neglecting (crossed sym- 
bols) and including (solid curve) the impurity-ladder vertex 
correction. The dashed curve is the line on which 04(0) = 0, 
H/H°^{Q). 
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planes, IT = — iV^ + 2|e|A, and the operator K2 in the 
pure limit (r^^ — > 0) is simply given by 



i^2(n,g,) = 2nN{0)T^Dd{2e), 

(kppGe„(p)G_,„(n + g3z-p))^ 



Dd{2e) 



ttN{0) 



(7) 



Here, the notation J = J l^A <f was used. Here- 



>p - J I27F 

after, a bracket {■ ■ ■)a with subscript a means an average 
over a, so that (• • ■)„ implies 2~^ J2a- ^^"^ impure case, 
there are also contributions from the impurity-ladder ver- 
tex corrections, and Dd{2e) should be replaced as follows: 



Dd{2e) ^ Z?d(2e + 1/T) 
+ Ed{2e) 



T 



l-TDs{2e) 



Ed{2e) 



2e->2e+l/T 



where F = {2t)~^ , and Ds and Ed are defined by re- 
placing lifpp in eq.(jZ|) into 1 and Wp respectively. In 
the dirty limit with 2'ktTco <C 1, the impurity-ladder 
vertex correction expressed by the second term of eq .(|S|) 
becomes important, although the SC phase itself in the 
present d-wave case is simultaneously suppressed. As is 
shown later, however, the MF-FOT appears only when 
2ttTcoT > 10, and hence, we focus here on this moder- 
ately clean case. In fact, as given in Fig^ a change in 
the i?c2 (T)-curve due to the inclusion of the impurity- 
ladder vertex corrections is extremely small even when 
(27rTcor)~^ = 0.1. Based on this result, the second term 
of eq.lIHJ and any contribution including the impurity- 
ladder vertex corrections in other terms of GL functional 
will be neglected hereafter in the text and in Appendices. 

We should note here that, strictly speaking, the eigen- 
states of the operator K2 in the dx2_y2-weLve pairing are 
not the LLs, and that there are nonvanishing off-diagonal 
matrix elements between LLL and higher LLs with in- 
dices of multiples of four. However, in the range of Maki 
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parameter considered in this paper, the instability Hne 
for the = 4 LL modes, defined by a4(0) = in our no- 
tation used below, lies far below Hc2{T) (see the dashed 
line in FigQ]) , and hence, we can neglect the off-diagonal 
elements in considering vortex states in = (and 1). 
Then, Wp in eq.Q may be replaced by 1, and our analysis 
using the LL basis becomes straightforward. When focus- 
ing on a projection Ag^-* onto the A^-th LL of Ag, (rj^), 
the corresponding eigenvalue of K2 is given by 



K2,N{qz) 



27tTN{0) dpf{p)CN(£^) 

X exp(-pV44)Jo(2Jsin(^)p), 



(9) 



where Cn is the N-th Laguerre function, j7o is the zeroth 
Bessel function, th = rnrn/kp, and the function / is 
defined by 



-p/t 



cos(2/p) 



sinh(27rrp) 



(10) 



The procedures leading to eq. © will be explained in Ap- 
pendix A. After, like in H = case, eliminating the high 
energy cut-off by defining Tco, we obtain the final expres- 
sion for the quadratic free energy 

T2 - A^(O) EE/ rf'^^«A'fe') |A(f Hr^)P, (11) 



N=0 9, 



where 



/•OO / 

aN{ql) = HT/T,o) + 2ttT J dpi{smh{2TTTp))~^ 

- .np)CN{^)e-(^y Jo{2Jsmi^)p)^. 

(12) 

For several purposes, it is convenient to expand aN{q'^) 
in powers of sm^{qzs/2): 

aN{ql) - aiv(O) + aj^^Q' + aj^^g^ + • • •, (13) 

where Q = Jsin((7^s/2). When the SC transition in the 
MF approximation is of first order and occurs within the 
iV-th LL, this MF-FOT line, which is a part of Hc2{T), 
lies in the region where aAr(O) > 0, and the quartic and 
sixth order GL terms have to be considered to deter- 
mine the Hc2{T)-\me. A possibility of an instability to 
an FFLO-like helical vortex solid in the A^-th LL may be 
studied, at least when |aAr(0)| <C 1 (see discussions below 



ea. (|45|l ). by focusing on 



and a|^\ where 



''N 



2TrT / dpp'f{p)e 
Jo 



1^^^ - 



As far as a^' > 0, the instability point to an FFLO solid 

state is determined by aj^'' = 0, and the corresponding 
transition is of second order. 



(14) 
(15) 



B. Quartic term 

The corresponding analysis for higher order (quartic 
and 6th-order) terms of the GL functional is more com- 
plicated than that for the quadratic one. As already ex- 
plained, the impurity-ladder vertex corrections will be 
neglected in the ensuing analysis. Hereafter, it is con- 
venient to work in a fixed Landau gauge A = (0, Hx, 0) 
and to represent the pair-field in the A'-th LL in terms 
of the corresponding LL orbitals UAr,/c(r^) 



Ar(r,) = 



1 



->H 



k 



(t>N,k.j UN,k{'^±), (16) 



where Sh = rnLy-K^^^. In the present gauge, Uiv,fc(rj_) 
is given by 



ujv,fc(r±) = 



\^(x+kr%f+iky 

+ -.e ^ " , (17) 



where we introduce the creation and annihilation opera- 
tors of LLs 



7r± = -|(n,±in,). 



(18) 



Hereafter, let us focus on a vortex solid within the LLL- 
subspace. The corresponding analysis in higher LLs will 
be given elsewhere. Then, the quartic term of the GL 
functional can be written as 







X 




X 


(a) 


(b) 


(c) 



FIG. 2: Diagrams expressing the quartic term of the GL 
functional. The solid line implies the Green's function G, and 
the impurity (dashed) line carries (27r7V(0)r)~^. 
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^4 = lY.Jd'rM{n^}){^f\r^lK\r^3)r^f{r^2)Af\r^,) 



1 v-^ f (Pr± - 

Z2 ^0M,J^0,k2,j(l^0M.j(l^0M.J J ^^^4({ni})uoj,^(r_Ll)MOfe2(ri2)MSfe3(''-L3)"Ofc4(r±4) 



(19) 



.(20) 



where 11^ = Il(ri). In the impure case K4 consists of three terms represented in Fig|5| and will be expressed as 
K4 — K^g, + Kihc- The term K4^^ is given by 



= 2ttTN{0)J2 



(2i£, + V • n*)(2i£, + V • n2)(2i£^ + V • n;) / ^ - 



(n2 ^ n4) 



» 3 3 

2TrTN{Q) / Y[dp,f{Y.P^ ^|^p|4gi(piv.nt+p.v.n.+p3v.n;)^^ ^ (jja ^ n4), 



(21) 
(22) 



where v = kpp/m, the function / is defined by ea. H10() . and the bracket ( )p implies the angle-average over the unit 
vector p = p/kp on Fermi surface. The sum of Fig. [21(b) and (c), K4bc, is given by 



4bc 



2t:T 



2t:T 



(2ie^ + V • n*)(2i£, + V • Hs) / - \ (2i£^ + v' • n*)(2ie, + v' • n4) / -,/^ 

+(n2 ^ n4) 



(23) 



4 4 

N{0) J f[dpJiJ2p){\wp\^e''^''^-^'^''^-^'^).{\wp'\^ 



i=l i=l 

where v' — kpp' /m. The following formulas which are derived in Appendix B are quite convenient. 



(25) 
(26) 



where A = pC* /'''h and C — Px + Wy is the complex coordinate. Using this identity and ea. l|()U|) in Appendix B, we 
obtain the following results 



X4a({n,})M*. (rii)uofe2 (r±2Xfc3 (i"-L3)uofc4 (rx4) 



Me'' 



X "ofci (r_Li)uo/c2 {^±2)uok^ (ri3)Mofe4 {^±4) ) 

' p 



+ (fc2 ^ fc4) 



27rrA^(0) 
V2 



3 3 

. /n^/'^/(E/') (4(ReC^)^/4({A,}) 

. — 1 A — 1 



+ C.C. 



(27) 



A'46c({n,})u5fe^ (r±i)wofc2 (i"i2)uSfc3 (ri3)wofc4 (i"±4) 



Sh 

2ttTN{0) 

T 



n2) 



"Ofei (r±i)woA;2 (ri2) 



i=l 1=1 

|2„i(p3v'-n3+p4v'-n4; 



Uok. (r±3)uofc4(i"J-4) 



p' 



+ (/C2 <-> fc4) 
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(a) (b) 
FIG. 3: Diagrams expressing the 6th order term of the GL functional. 



2TrTN{0) 
tV2 



„ 4 4 

^ . — 1 — 1 



i=l i=l 



Al,2=Pl,2C*/'^ff; -'>3,4=P3,4?*/th ' P,P' 



+ C.C. 



(28) 



where the suffix \ij = PijC /'''h impUes = piC /th and = pjC,* /th- Here the function I-^ for the quartic term 
is given by 



,4 

In (/4({AJ)) - -T E l^'l' - «(^l3 + AL) --AK+ A3)(A2 + A4) - ^(fci3Al3 - fc24A24), 



(29) 



where S, — p'x ~^ '^p'y 1 kij = ki — kj , and Xij = Xi — Xj . Finally we have a quartic term 

= j^i^} J2 J2 '5fcl+fe3,fc2+fc4"t^4({fc»})e-^('=" + '='^VS.fcl.J</'0,fc2.j'/'S,fc3 j'/'0.fc4 J, 



V4{h}) - 2^r / l[dpj{j:LiP^) ■ (4(ReC')*/4({A,}) 

/ nrfP^/(i:tipO • (4(RcC^)^(Re^^)^/4({AJ) 



(30) 



A4=0; Ai=piC*/Tff / p 



Al,2=Pl,2C*/'''ff ; A3,4=P3,4^* /rff 



) + C.C. (31) 
' p,p' 



We will show later (see Fig. 5 (a)) that, consistently with the neglect of the second term of eq.®, the second contri- 
bution to eq.|j2U arising from i^4bc is safely negligible compared with the first K^^, term in the relatively clean case 
with 27rTcoT > 10 of our interest. 

C. 6th order term 

When we restrict the pair-field into the LLL subspace, the 6th order term of the GL functional are expressed as 
follows 



^6 = -iE/^'-^^6({n.})(Af (r^OAf (rx3)Af (r^5))*Af (r^2)Af (r^4)Af (rx6) 



(32) 



002 Z-^ •rv,Ki,j-rO,k2,j(Po,k3,jVO,ki,j(Po,k5,j'PO,ke,j 



X / -g^Ke{{U.i\)ulf^^ (r±i)Mofe2 (r±2)w5fc3 (ri3)Mofe4 (r±4)wSfc5 (i"-L3)wofc6 (i-±4) 



(33) 



In contrast to the quartic term, the kernel Kq also includes diagrams (see Fig|31(b)) with two or three impurity lines 
in addition to those with a single impurity line such as Fig[21(b). Fortunately, according to the statement following 
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FIG. 4: Results on #4 = 2J"4/[iVij^y^^r2iV(0)r^({|A'°' p))^] calculated under three different conditions and in the pure (tTco = 
00) limit. In the open-square symbols, the square vortex lattice is assumed with the nonlocal contribution in F4 included, 
while the triangular lattice is assumed in both the crossed symbols with the nonlocal contribution and the solid curve with no 
nonlocal one. Note that T* at which _F4 = is insensitive to the nonlocal contribution and to the types of lattices. 



ea. (|31|l . all terms other than Fig|21(a) may be neglected in the range of purity parameter we focus on. The diagram 
FigEl (a) is expressed as 



Kg = 27rTA^(0)^ise 



Z6 + + Z2 + 



Z3Z6 



Z1Z4 



Z2Z5 



Z2- Z3 + Z4 Z3 - Z4 + Z5 Zi - Z2 + Z3 



(34) 



where Zi = 2i£o- + v • 11^ for even i and 2iia- +v • 11* for odd i. It is easily seen that, due to the symmetry with respect 
to Hi and 11*, the above expression can be represented in terms only of, e.g., its first and fourth terms. 

Firstly, let us calculate the contribution to the GL fimctional of K^a- Using the parametric representation (see 
Appendix A), it is written as 



i4:6a({n,})M*fc^(ril)u0fc2(l"±2)MSfc3(l"-L3)M0fc4(r±4)uSfc5(rj.3)u0fc6(l"-L4) 

5 5 



(35) 



r I i— ►r 



-27rTiV(0) J l[dpJ{Y,P) J 



1=1 i=l 



27rTiV(0) 
V3 



5k 



ki+k3 + ks,k2+k4, + ke 



Wo,fei(l"il)wo,fe2(r_L2)Uafc3(l"i3)wo,fc4(l-J_4)Uo,fc5 (l"±5)wo,fc6(l"±6) 
5 5 

f[dp,fiY^p) (8(ReC^)6/6({A,}) 



i=l 



>'6=0, \i^G=PiC* /th I P 



(36) 
(37) 



where the function Jg is given by 



In 



2—1 z:odd i:cvcn z:odd i:cvcn 

(-i j'):odd (i j'):cvcn (i j'):odd ('ij'):cvcn 



(38) 



Next we examine K^d- Using the parameter representation, this term is expressed by changing the integral variables 
in the above expression as /32 — >■ P2 + Pa, Pa — ^ —pz-, and P4 — > pa + P4. 

In this way, we can write the 6th order term of the GL free energy functional in the form 



^6 



^(0) 
3^35*1 



E E '^'=i+'=3+'==''=2+fc4+fc6''4>({fcj}) e ^< 



''0,fci,iV^0,fc3j''^0,A;5,jV'0,fc2,i'P0,fe4j'P0,fe6j' 



(39) 
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(b) 

FIG. 5: Numerical results of the dimensionless coefficients (a) V4 — V4{{kj — 0})/r^ and (b) Ve = Ve{{kj = 0})/r|f on the 
_H"c2-curve at lower temperatures T < O.STco- The value {2-kTcot)~^ — 0.05 was commonly used. In (a), the cross symbols 
represent the result due only to Fig|2ja). 

where the summation J2(i m) taken over the pairs {(1, 3), (3, 5), (5, 1), (2, 4), (4, 6), (6, 2)}, and Vq is given by 
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FIG. 6: High H and low T region of MF H-T phase diagrams in (a) clean limit ((27rTcoT) ^ = 0) and (b) the moderately 
clean case ((27rrcoT)^^ — 0.05). See the text for further details. 



In deriving the MF phase diagram and its impurity 
dependence, we will use an additional approximation be- 
low. Since a nonzero magnetic field plays the roles of 
cutting off the low T divergences of coefficients of the 
higher order GL terms, the orbital depairing effect aris- 
ing from the gauge invariant gradients Tlj has been in- 
corporated nonperturbatively in the above expressions. 
Instead, algebraic kij dependences have arisen in the ver- 
tices Vm{{kij}) with m = 4 and 6. On the other hand, 
the Gaussian fc^ -dependences in J^^ and JFg are direct 
consequences of restricting the pair-field into the LLL 



subspace and also appear in the familiar GL expression 
with spatially nonlocal higher order terms. That is, the 
additional kij dependences in Vmi{kij}) can be seen as 
spatially nonlocal contributions to the higher order GL 
terms and affect the structure of vortex solid. Actually, 
in LLL and the case with no paramagnetic depairing, this 
nonlocality in the quartic GL term results in the struc- 
tural transition between the rhombic and square vor- 
tex lattices 20j. However, an energy difference affecting 
the lattice structure is extremely small reflecting a weak 
structure dependence of the Abrikosov factors (denoted 



9 



as (3a and 7^ below). For instance, as shown in Fig^ 
the nonlocal correction to F4 and thus, to T* is negli- 
gibly small. Therefore, at least as far as the SC transi- 
tion in the MF approximation at is concerned, such 
nonlocal corrections are safely negligible. For this rea- 
son, the local approximation for the higher order terms 
will be used hereafter, and Vm{{kij}) will be replaced by 
Vm = Vm{{kij = 0}). Then, the GL functional derived 
microscopically takes the form 



^loc = iV(0) / 



Temperature variations of the coefhcients V4 and Vg 
calculated along the 77c2 (T')-curve are shown, respec- 
tively, in Fig[Sl (a) and (b). To clarify that the contri- 
butions of Fig[21 (b) and (c) are safely negligible, V4 in 
FigEl (a) was calculated in terms of {2TrTcoT)~^ — 0.05, 
the value used in FigEl (b) below. The coefficient V4 is 
negative at lower temperatures, while Vg is positive over 
a broad region so that the GL expression H41|l with a 
MF-FOT at Hc2{T) is well-defined. 



D. Mean field Phase diagram 

Below, the MF phase diagram will be examined based 
on the functional, ea. (|lT|l . First, let us neglect a possi- 
bility of an FFLO-like state and assume a straight vor- 
tex solid independent of j as the MF solution. Then, 
the first term in the bracket of ea. H41() is replaced by 
X]j '^0(0)1 Aj p, and the MF solution is obtained in a stan- 
dard way. The character of the MF transition at Hc2{T)- 
line changes with increasing field from the second order 
one to a discontinuous one at a temperature T* where V4 
becomes negative (see Fig0]and Fig|5l(a)), reflecting that 
the spin depairing is more effective upon cooling and with 
increasing H. To obtain MF results in T < T*, higher 
order terms are necessary in the GL expression. Accord- 
ing to Fig[31(b), the coefRcient Vg is positive over a broad 
temperature range, and thus, the expression ea. (|41|l ter- 
minated at the 6th order term will be sufficient for our 
purpose. Further, let us introduce the effective coeffi- 
cients b = V4PA of the quartic term and c = VqJa of the 
6th order term, respectively, where 



|Af (rx)l^ 
((|Af 

|Af (rx)l^' 

(M¥y) 



(42) 




0.2 0.3 

T/Tco 



. (41) FIG. 7: Temperature dependences of c^Q^ and a^' on Hc2{T) 



-(2) 



and near Tfplo, where 5q = %" /T[f (s = 1, 2), and the pure 
limit was assumed. The corresponding curve of a^^^ = a^^^ /rfj 
is also shown for a convenience. 



Then, the following MF results in T < T* are found. 
First, the MF transition point in T < T* and in LLL is 
determined by the condition 



ao(0) 



16c 



(43) 



while the supercooling (superheating) point is given by 
ao(0) = (ao(0) = 6^ /(4c)). Next, the energy barrier 
^^barr between the | A*^°) | = solution and the jump value 

of |A("'| at the transition, |A(°)|c = y^3|6|/4c, is given 
by 



\b\^ 

t/.a.^iV(0)i^. 



(44) 



Further, by calculating the mean squared amplitude 
(|(5A('^)p) of the Gaussian fluctuation ^A(°) when 
ao(0) = ao,c and in 2D limit, one also finds 



|<5A(°)|' 



rc2|AW| 

^barr 



(45) 



where Tc2 is the MF transition temperature. Thus, the 
fluctuation strength is enhanced with decreasing |6| and 
increasing c and, as expected, is measured at Tc2 by the 
inverse of the energy barrier. Hence, if this MF-FOT 
occurs as a true FOT in real systems, a clear hysteresis 
is expected in a system with weaker fluctuation. 

However, in higher fields and lower temperatures where 
the spin depairing becomes more important, an FFLO- 
like helical vortex solid within LLL may become more 
favorable. As far as the width ao.c is sufficiently small, 
the Qz-dependent terms have only to be incorporated in 
the quadratic terms in A^'^''. That is, this structural tran- 
sition line Hfflo{T) between the FFLO-like solid and 
the straight vortex solid may be discussed within the co- 
efficient ao{ql). Actually, according to the calculation 
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results of V4 and in FigO ao,c in Fig|B| is at most of 
the order of 10~^. Assuming a helical state with the j- 



dependence A 



(0) 



gi9^js in _ff > i/pFLOj a second order 



structural transition line -ffpFLoCT^) is obtained within 
LLL according to ca. l(T^ as the curve a^\T) = if 
a^\T) > there. Numerical results on a^^ and a'^^ are 
shown in Fig[7| We find that a^^' near Hc2 is always pos- 
itive along the HpFi,o{T)-\me and increases upon cool- 
ing at a fixed field above Hfflo{T). Further, since the 
paramagnetic depairing effect is enhanced with increas- 
ing field and decreasing temperature, as the example in 
FigEl shows, a possible iJpFLO (^)-curve should decrease 
upon cooling. 

Now, the MF phase diagrams Fig |H| following from the 
GL coefficients derived above and for two different impu- 
rity strengths ((rTco)^ ^-values) will be explained. The 
used value of paramagnetic parameter /io-ff°2'^(0)/27rrco, 
corresponding to the Maki parameter, is 0.8, which leads 
to the value T* = OMTco in the pure hmit (Fig.6 (a)), 
where H°^^{0) = 0.56(/)o/(2<g) is the 2D orbital limit- 
ing field at T = in pure limit, and ^0 = kF/2'!TmTco is 
the coherence length. The field values in the figures were 
normalized by H°^^{0) (i.e., h = H/H°^{Q)). The curve 
H*{T) indicated by the cross symbols is defined by the 
condition V4 = and may not be directly seen in experi- 
ments. In contrast, the portion (open circles) in T < T* 
of Hc2{T) on which the MF SC transition is discontinu- 
ous is experimentally measurable together with (if any) 
the second order transition line i?FFLo(T') (solid curve) 
to the FFLO vortex solid. In the temperature regions 
where the MF-FOT does not occur, the higher one of the 
dashed curves indicated by = or = 1 becomes the 
i?c2(r)-line, on which ao(0) or ai(0) = 0, and the MF 
transition there is of second order. 

It will be important to, in relation to real phase 
diagrams of related materials, understand how the 
-ffpFLO (T) curve and the characteristic temperatures are 
affected by the impurity strength. By comparing Figs. 6 
(a) and (b) with each other, the region ifpFLoC?") < H < 
Hc2{T) is found to be easily lost by a sHght increase of 
impurity strength (tTcq)^^. In contrast, the onset T* of 
the MF-FOT behavior is relatively insensitive to the sam- 
ple purity. Nevertheless, when {2ttTcot)~^ goes beyond 
0.095 while the value hoH°^{Q)/2ttTco = 0.8 was kept, 
the MF-FOT region at Hc2{T) is also lost, and the MF 
transition at Hc2{T) is continuous at all temperatures. 
This result is contrast to other works 0, in which 
the presence of a MF-FOT was argued under the use of 
dirty limit. We find that, instead, the FOT obtained in 
the dirty limit |l2l | never occurs in T — )■ limit when 
EpT > 1 under which the usual dirty limit may be valid. 
On the other hand, the results in ref.0 are derived by 
completely neglecting the orbital depairing and are not 
comparable with the present ones. Further, we stress 
that, in contrast to results given in a previous work [l]| 



taking account of both the orbital and spin depairing ef- 
fects, the results in Fig|Bl imply that always T* > Tfflo 
in the present quasi 2D case under a perpendicular mag- 
netic field. 

Since, as already mentioned, the width ao,c in FigElis 
relatively small, the MF-FOT there may be regarded as 
being weak. However, it does not mean a strong fluctua- 
tion. Actually, in systems with a large A^(0) in zero field 
such as CeCoIns, the fluctuation strength T/Uhair itself 
becomes extremely small in the low T region of our inter- 
est. On the other hand, the magnetization jump value 
AMc at the MF-FOT should be quite small compared 
with the applied field in order to justify our neglect of 
a spatially varying internal magnetic field. In CeCoIn5 
under an applied field in the tesla range, this condition 
is well satisfied Hi (see also Sec.IV). 

In T < Tnoxt where ai(0) < ao(0) (< 02(0)), the 
Hc2{T) line and hence, the vortex lattice itself just be- 
low it are determined by the next lowest (A^ = 1) LL. 
According to the a^^^(T)-curve in Figd this A^ = 1 state 
is not a modulated FFLO-like state but a straight vor- 
tex solid. Further, note that Tncxt and Tfflo are close 
to each other (see Fig|Hl(a)). These results imply that, 
in quasi 2D materials under flelds perpendicular to the 
SC layers, the presence itself of the FFLO-like vortex 
solid is very subtle even in pure limit. Thus, a competi- 
tion between the FFLO-like solid within LLL and a solid 
within the A^ = 1 LL has to be examined just above the 
HpFhoiT) line. Since this is an issue of a transition be- 
tween vortex lattice structures defined within the planes 
perpendicular to the field, a detailed description of the 
stable vortex lattices in d-wave pairing cases is required 
to address this. As already mentioned, however, the non- 
locality, at least, of the sixth order GL term affecting the 
in-place lattice structure was neglected in this paper. We 
will postpone a study of structural transitions to higher 
LL solids in H > i?FFLO to our future study. 



E. Properties of Helical Vortex Solid 

Here, we briefly comment on linear responses and 
elastic properties in vortex solid phases, primarily in 
an FFLO-like helical solid with a modulation Aj ~ 
Ag^e'*"^* in LLL with g.,„ 7^ 0. To examine the electro- 
magnetic linear responses in an ordered vortex lattice 
phase, we have only to focus on the gradient terms with 
an external gauge fluctuation substituted and to examine 
the Gaussian fluctuation around a MF solution of Aj . An 
appropriate form, consistent with the above microscopic 
analysis, of the gradient energy will be 



rad 
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A(ri,z), 



(46) 



where a continuous variable z was used for the coordinate 
parallel to the field. If considering the SC response in the 
X or y direction, 11 needs to be accompanied by a gauge 
fluctuation a = + in the form 11 — aj^. Within 
the Gaussian approximation for the fluctuation, no cross 
terms like aj_ Uz appear because any term cx aj_ becomes 
off-diagonal with respect to the LLs and hence, zero af- 
ter spatial averaging. Hence, the linear responses in the 
parallel and perpendicular directions can be considered 
independently. 

First, let us consider the response parallel to the field 
in which = 0. If the spatial variations perpendicular 
to H of the MF pair-field solution are described within 
LLL, the argument 11^ in A and B can be replaced by 
r^^. For the moment, let us focus on the helical solid 
phase in which B{r~^) > 0. Assuming the fluctuation of 
Aj to be dominated by that of its phase 0, the fluctuation 
part of 5.?"grad simply becomes 



-2dz^az + {dl 



(47) 



Here, 5<f) is the phase fluctuation, and higher gradient 
terms were omitted in the expression following the sec- 
ond equality. Further, the fact that the MF value 
of is Bir'^) / A{r~^) was used. As seen in the semi- 
nal work the same form as above occurs in the case 
with no orbital depairing in an isotropic 3D system (see 
eq.(31) in ref.j^). This familiar result (|47l) implies that, 
in the helical solid state, there is no stationary current 
even along H and that a supercurrent cx can flow in 
this direction, although the superfluid rigidity cx B{r^) 
vanishes with approaching i?FFLo(2^)- 

The absence of SC response in x and y directions can 
be seen, as well as in the ordinary (g^ = 0) case, most 
easily within the phase-only approximation j2^. Not- 
ing that any term accompanied by the in-plane periodic 
variation in a vortex lattice disappears after averaging 
spatially and that the uniform fluctuation part of d±S(j) 
is given by r^^ (z x s) , one finds that the fluctuation of 11^ 
in the limit of vanishing wave number can be expressed 
everywhere in A and B as r^'*((z x s) — aj^)^, where s is a 
uniform in-plane displacement of a vortex lattice. After 
integrating over s, aj_ disappears together in the resulting 
fluctuation free energy, implying no SC response. This 



non-SC perpendicular response is a consequence of the 
establishment of the Josephson relation aj^ = (z x s) 
(leading to an electric field E proportional to a vortex 
velocity v^, i.e., E = — x H) and is quite different 
from the corresponding one in the case with no orbital 
depairing in the ideal isotropic 3D system which is a con- 
sequence of spontaneous formation of a direction of the 
helical modulation jl^ . 

Next, let us briefly discuss the vortex elastic energy. 
The shear elastic energy for in-plane shear distortions is 
obtained, as in the case of straight vortex lattice, from 
the phase fluctuation energy. It was already shown in 
ref.|22| that the quartic dispersion on the 2D wave vector 
qj_ of the massless phase fluctuation and hence, the way 
of identifying the shear elastic mode with the phase fluc- 
tuation are unaffected by any change of the higher order 
terms of the GL model. No calculation of shear modulus 
will not be performed here, and we simply assume, just 
for qualitative consideration, the resulting shear energy 
term to be positive definite. In contrast, some comment 
will be necessary on the tilt energy which should change 
through Hfflo(T). For simplicity, we focus on its ex- 
pression in type II limit with no internal gauge fluctua- 
tion 0,1131. According to ea. l|T7jl . the tilt energy in type 
II limit takes the form 



E< 



tilt 
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where \kz\ = \\qz\ — qm\, and the relation E3.^2^ 
5(t>q — -~i{q]_^r^'^){qj_ x s^)^ between the transverse dis- 
placement S"'" and the phase fluctuation was used. On 
the other hand, in the straight vortex solid (i.e., in H < 
HpFhoiT)) where B{r^) < 0, the dispersion k1{k1 + 
4g^) in eq.gHIl is replaced by 2\B{r-^^)\ql/ A{r-^^). Thus, 
the macroscopic properties in both the vortex solids are 
qualitatively the same as each other. 

However, the tilt rigidity decreases with approaching 
H-pFwiT). In particular, the k1/q\ dependence of the 
tilt modulus just on HjfphoiT) implies a s/iort-ranged 
phase coherence of the vortex solid there which is consis- 
tent with the vanishing superfluid rigidity on i/pFLO (T) . 
Because this absence of SC order in the critical vortex 
solid (i.e., the case just on i?FFLo(T')) is due to a soft- 
ening (or weakening) of an elastic modulus, an enhanced 
peak effect due to an increase of the critical current is 
expected near HpFi^o(T), as well as the ordinary one 
near Hc2(T) |27l |. in real materials with random pinning 
effects. This conjecture may be useful for searching a 
transition curve to the FFLO-like vortex solid. 



III. STUDY OF GENUINE PHASE DIAGRAM 

In this section, the real phase diagram of systems de- 
scribed by ea. l41|) is studied by including the fluctua- 
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tion effects. We will focus here on the range H*{T) < 
H < i?FFLo(r) and hence, rewrite the first term in the 
bracket of eq.^ into E,(«o|Af ^ + 7o|Af - A^^J^) 
with 7o > 0. Although we have not extended our sim- 
ulation work to (if any) the range H > i?FFLO defined 
within the lowest LL, fluctuation effects similar to those 
in if < -Sfflo should be also expected in such higher 
fields (see below). The partition function we should ex- 
amine is 

Z = Tr^cxp(-J^), (49) 
where the functional J-' = J-ioc/T is rewritten as 



(50)' 

where ^E* j (r) is the rescaled order parameter field defined 
within LLL, 7 > 0, 



r|7V(0)^ 
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1/3 
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(51) 



and the in-plane coordinates r were normalized like 
J^/^H r. Further, Hq denotes iJc2(0) in the case with 
no MF-FOT. Since, as mentioned earlier, ao,c measuring 
the difference {Hc2—Hq)/Hq is small for the Maki param- 
eter value used in this paper, Hq will not be distinguished 
from H(.2 below. Note that, except a numerical factor, 
is identical with Uhm-r/T. The T-dependent param- 
eter measures the fluctuation strength, while the 
magnetic field dependence is assumed only in a. Thus, a 
change of can be regarded near Hc2 as a change only 
of T in eq. H49|l under the same J-ioc with fixed values of 
ooi 70: Va, and Vg, if the magnetic field variable a is ap- 
propriately rescaled, and a difference in the anisotropy 
(i.e., 7o-value) plays no important role. 



A. Theoretical consideration 

First, let us explain how an ordinary physical picture 
on low energy fluctuations and their effect on orderings 
lead to the absence of a g enuine transition at Hc2- Like 
in the familiar case [sl, |21| with a second order MF tran- 
sition at Hc2, let us start from a description deep in the 
ordered phase. First, since an inclusion of the orbital de- 
pairing requires the presence of field-induced vortices, a 
low energy excitation in the ordered phase is inevitably 
an elastic mode of a vortex solid. As already noted in 
the final part of Sec. II, it is clear that the form of the 
elastic energy and the relation between the phase fluc- 
tuation and the vortex displacement |2^ are essentially 
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FIG. 8: Two candidates of a schematic H-T phase diagram 
of bulk type II superconductors with paramagnetic depair- 
ing in the case with thermal fluctuation and with no vortex 
pinning effect. The solid curve, the thick dashed curve, and 
thin dashed one are the melting line Hm{T) of vortex solids, 
the MF-FOT curve on Hc2{T), and the ordinary Hc2 curve of 
the second order MF transition, respectively. Dashed curves 
are not genuine transition lines. A possible presence of the 
structural transition to an FFLO-like vortex solid, neglected 
here for simplicity, does not lead to any essential change on 
Hm{T) and H,2{T) 



unaffected by differences in the forms of higher order GL 
terms. Then, it is clear from the previous works |2ll |23| 
that the phase fluctuation is marginally relevant even in 
3D case, and that the rigidity controlling the quasi long- 
ranged phase coherence in the vortex solid is the shear 
modulus. '2l'| That is, if the vortex solid is melted at Hm 
below Hc2 , wc have only short-ranged orders for both the 
phase and the vortex position in the vortex liquid above 
Hm with no flnite shear modulus, and thus, the vortex 
liquid should be continuously connected with the normal 
phase above Hc2 In this sense, the MF-FOT at 

Hc2 should not occur. Further, to understand this from 
another point of view, let us note that the quasi 2D SC 
order parameter in the lowest LL has the form 



^{tz)^A{z)e-y n (C-e.(^)), 



(52) 



4=0 



where = x + iy, z = js, ^i{z) is the complex coordi- 
nate perpendicular to H of the i-th vortex, and a Landau 
gauge was assumed for the external gauge field. Since 
the vortex positions are highly disordered above Hm, the 
fluctuation effect above H^ is essentially described only 
by the amplitude A{z). However, A{z) depends only on 
z irrespective of the form of the higher order GL terms. 
That is, since the amplitude fluctuation itself has a re- 
duced dimensionality and is ID-like in 3D systems Eoj 
even in the present case, the amplitude fluctuation will 
push the vortex lattice freezing field H^ down to a lower 
field, and the MF-FOT should not be realized as a true 
FOT in real 3D systems jSa]. 

Once noting that the g^-form of the tilt elastic term in 
H < i/FFLO is replaced in H > -ffFFLO by (Igr^l - g^)^, 
the above argument precluding a genuine FOT at Hc2 



13 



can be applied to such higher fields with no modification. 
Note also that the continuous transition at TJfflo is a 
structural transition between the ordered vortex solids 
and hence, together with FFLO-like states themselves, 
will not be lost due to the SC fluctuation as far as H„i > 
HpFLO (see Sec. IV). 

The theoretical discussion given above implies that, as 
far as Hm < Hc2, a genuine transition at Hc2 cannot 
occur. Then, it is at least natural to expect that, fol- 
lowing the familiar case in T > T* with a second order 
ifc2-transition, < Hc2 at any nonzero T and hence 
that, as in Fig|Hl(a), a genuine transition at Hc2 never oc- 
curs at any nonzero temperature. However, caution will 
be necessary in the present case with a MF-FOT. For 
instance, a Hm defined in terms of the Lindemann crite- 
rion may lie above Hc2 in a case with weak enough 
fluctuation because any elastic modulus, proportional in 
the lowest LL to some power of (I^Pp), is nonvanishing 
on approaching Hc2 from below. Of course, there will be 
no possibility that the actual Hm lies above Hc2- Then, 
one may consider another possible phase diagram FiglHl 
(b) realizable for a case with weak enough fluctuation, in 
which, reflecting a reduction of thermal fluctuation upon 
cooling, HmiT) at low enough temperatures virtually co- 
incides with Hc2, and both a large jump of magnetization 
(reflecting a large condensation energy) and a tiny hys- 
teresis, arising from the vortex lattice melting, will be 
seen at Hc2 on the solid curve. An important point is 
that, even in the scenario Fig|Sl (b), a large hysteresis 
occurring as a consequence of a large jump of magneti- 
zation, i.e., a hysteresis at thermal equilibrium resulting 
directly from the MF-FOT, cannot occur. According to 
the experimental data [13, lla ll^ in the situation of our 
interest in this paper corresponding to CeCoIns under a 
field H II c, the phase diagram of Fig|S| (a) seems to be 
always realized. However, it is theoretically valuable to 
clarify whether FiglHl (b) may be realized or not in real 
systems. Results of our numerical simulation, performed 
for examining our theoretical arguments given above, will 
be reported below. 

B. Simulation results 

In this subsection, we explain our methods and results 
of Monte Carlo (MC) simulations for the model ea. (|49|l . 
Our simulation method closely follows that used in the 
simulations [3^l8^ for the familiar GL model with a pos- 
itive quartic term in place of — 1/3| in ea. H5()|l . On a fixed 
SC layer, the pair-field ^E" is expanded in terms of the 
LLL basis functions (pnjx^y) consistent with a quasi pe- 
riodic boundary condition |32| as 5" = '^n'^n4'n{x,y), 
and the system sizes and Ly of a rectangular cell 
satisfy the commensurability with a triangular lattice 
ground state through the relation L^/Ly = ^/3Nx/2Ny 
(= 2Nx / VSNy) in 2D (layered) case. Note that, as men- 



tioned earlier, the ground state in the present case is, 
due to the neglect of nonlocality in the GL higher order 
terms, a triangular lattice, although the pairing state as- 
sumed originally is a four-fold d-wave one (see the sen- 
tences above eq.(0TJ). In the layered periodic 
boundary condition is assumed across the layers. The 
system sizes we have studied were {Nx, Ny) = (6, 4) and 
(8, 6) in 2D case and (6, 6) in the layered case. The 
Markov chains for c„ are generated by the Metropolis al- 
gorithm. For 2D case (layered case), we used 5 x 10^ 
(1.5 X 10^) MC steps for thermalization and another 
5 X 10^ MC steps for further observation in both cases. 
As in the figures in Sec. II, we have assumed the value 0.8 
for ^o^fc2''(0)/27i"Tco. Further^hen assuming a N{0) 
value appropriate to CeCoIus jsj and that T/T^q = 0.1, 
the |/3(T)| value is estimated to be in the range between 
2.0 and 3.0 used in the simulations. 

To study fluctuation effects on the MF-FOT, the mean- 
squared average of the pair-fleld (|5'P) was calculated. It 
corresponds to the magnetization when ao is the measure, 
primarily, of H . Hence, if it shows not a true discontinu- 
ous jump but a rounded behavior at Hc2 broadening with 
decreasing |/?|, a genuine FOT at Hc2 is judged to be 
absent. Further, as a measure of the vortex-positional 
ordering (vortex-solidiflcation), we have examined the 
structure factor ^(k) deflned [23 as the Fourier trans- 
form of the correlation function of |^'(r)p. In the figures, 
a-dependences of these physical quantities are shown. 
Although a change of a implies a change of H at a fixed 
T, the vortex density is fixed at any simulation. 

First, let us present and explain 2D simulation results 
in which 7 = 0. The obtained results (symbols) are 
shown in Figr^and llOl The corresponding MF curves are 
also drawn for comparison. The feature, that the simu- 
lation data lie above the MF curve in Hm < H < Hc2, is 
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FIG. 9: Numerical data of a-dependence (i.e., H- 
dependence) of (j^'P) in 2D case, for |/3| = 2.0 (crossed sym- 
bols) and \(3\ = 3.0 (plus symbols). Dashed and solid curves 
are the corresponding MF results. According to eg. 14311 . the 
MF transition point corresponds to q = 0.75 (1.69) in (3 — 2.0 
(3.0). The system size (6,4) was commonly used. 
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(a) 




FIG. 10: S(k) data corresponding to the [/?[ = 2.0 data in 
FigHat (a) a = 0.75 and (b) -2.0. (c) S'(k) at a = -2.0 for 
larger system size (8, 6). 

not surprising but presumably a reflection of a dimension- 
ality dependence of the amplitude fluctuation (compare 
with Figll II below). It is found in the literature js^l that 
a similar dimensionality dependence appears in the mag- 
netization in the case with a positive quartic coefficient. 
As is clear particularly from the |/3| = 2 data of Fig|51 
the nearly discontinuous jump of (I^Pp) at Hc2 is rounded 
due to the fluctuation, and thus, no genuine FOT has oc- 
curred at Hc2. We note that ao{T = O.lTco) ~ 50 for the 
parameter values used here. Hence, if the abscissa in 
FigElis reexpressed as the reduced field [H — Ho)/Ho, 
even this rounded behavior of (l^'P) for \(3\ = 2 cannot 
be distinguished from a strictly sharp discontinuity, and 
the difference Hc2 — Hm will not be visible (see Sec. IV). 
Thus, the presence of a genuine FOT cannot be argued 
through merely a steep growth of in a real system 

with weak SC fluctuation. 

Next, let us examine whether the melting position co- 
incides or not with the MF transition field Hc2{T). The 
melting transition is widely believed to be a weak FOT, 
and this should be found [22| in Monte Carlo simulations 
as a tiny discontinuity in thermodynamic quantities. Un- 
fortunately, due primarily to numerical difficulties, our 



simulation is restricted to too small systems to observe 
such a discontinuity. For our purpose of determining Hm, 
however, it is sufficient to find where the Bragg peaks of 
the vortex lattice disappear. Figure [T7)l shows snapshots 
of the structure factor on |^(r_L)p when |/3| = 2.0. No 
vortex positional ordering is seen at Hc2- By comparing 
FigEl (a) and (b) with Fig^l one finds that nearly sharp 
Bragg peaks appear at Hm below Hc2, while most of the 
entropy has been lost near Hc2 above it. The two field (or 
temperature) scales, one characterizing the steep growth 
of (l^'P) and another corresponding to the sudden growth 
of vortex positional ordering, are clearly distinguished. 

We have also examined the size dependence of S{k.) 
data. By comparing the data in FigEII (b) and (c) with 
each other, it will be clear that the six- fold symmetry 
of Bragg peaks is more remarkable in the former, i.e., 
in a smaller size. It means that the solidification is en- 
hanced by the boundary condition in a smaller system. 
This size dependence will be sufficient for justifying our 
expectation that the vortex solidification (or melting) oc- 
curs below Hc2- 

The corresponding results in a layered (quasi 2D) sys- 
tem consisting of four layers are shown in Figslllland ll2l 
where the parameter values 7 = 0.25 and \(3\ = 2.0 or 3.0 
were assumed. The obtained computation results are es- 
sentially the same as in 2D case, except the feature that 
the difference Hc2 — Hm became narrower in the layered 
case. Thus, the anisotropy or dimensionality, i.e., the 
magnitude of 70, does not seem to induce an essential 
change of the true phase diagram. The nearly discontin- 
uous behavior at Hc2 is also smeared out in this quasi 2D 
case as the fluctuation is enhanced, and no hysteresis is 
present accompanying this behavior at Hc2 ■ The simula- 
tion results given above imply that, at least for |/?| < 3.0, 
the correct phase diagram is FiglHl(a). 

Next, we report on consequences of an extension of 
simulation for the layered system composed of four lay- 
ers to weaker fluctuation cases with |/3| > 3.0 (i.e., at 
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FIG. 11: Numerical data of (j^'l^) in the layered case com- 
posed of four layers. The crossed symbols (plus symbols) are 
results in \f3\ = 2.0 (3.0). 
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FIG. 12; Structure factor defined from the correlation func- 
tion of |*(r)p for = 3.0 at (a)Q = 1.65,(b) a = 1.25. 



lower temperatures). As the numerical data in Figll3l 
(a) show, a hysteresis in the vicinity of Hc2 suggestive 
of a genuine FOT appears between two curves for 

\P\ = 3.5, respectively, in increasing H (corresponding to 
a heating) and decreasing H (corresponding to a cooling) . 
However, this hysteresis is not due to the vortex lattice 
freezing or melting included in the scenario of Fig|Sl (b) 
because, as the simulation results in < 3 have shown, 
the hysteresis accompanying the melting is unobservably 
small (Although the melting field H„i is estimated to lie 
close to a = 2.0 through SCk) data, it is not easy to con- 
clude a separation between Hm and Hc2 for this case). 
Further, since, as mentioned below ea. 1)51(1 . an increase 
of |/3(r)| can be regarded as a decrease only of T in the 
partition function H49|) under the same J-'ioc (|41|) , it is dif- 
ficult to imagine a scenario that the MF result would 
become exact in a very low but finite T region. Actually, 
the hysteresis in Fig^l (a) is not due to a genuine FOT 
in thermal equilibrium: The a = 2.19 data in Figll3l (h) 
show that the system, at least in the vicinity of H(.2 , has 
not reached the thermal equilibrium even during the MC 
steps we could observe. 

Actually, a similar situation occurs, if the fluctuation 
is weak enough, in other systems with a MF-FOT but 
no true phase transition. To illustrate this statement, 
we show in Fig^J results, corresponding to FigEl for a 
familiar ID GL model with a negative quartic coefhcient 



= / dx[a\^{x)\'' + J 



dx 



(53) 



where is a function only of a;, and a in this case may 
be regarded as a temperature variable. This ID model 



3.0 



2.0 



1.0 





H-increase ■ 


H-decrease i 





1.0 



2.0 

a 

(a) 



3.0 



4.0 



S 1 




Monte Carlo step 
(b) 

FIG. 13; (a) Numerical data, similar to FiglTT] for |/3| = 3.5. 
(b) The history of (|*p) at a = 2.0 and a = 2.19 in the H- 
increase process. Note the strong dependence on Monte Carlo 
steps of (1^1^) when a — 2.19, i.e., in the vicinity of Hc2- 



is used here for comparison because the fluctuation in 
3D GL model within LLL is expected to be similar to 
that of the corresponding ID GL model in zero field. We 
chosen the values \'w\ — 5.0 and / — 4.0, together with 
the system size composed of 10^-sites. Due to a smaller 
number of degrees of freedom in SC fluctuations in this 
ID case, a tendency of a full relaxation to the thermal 
equilibrium (i.e., a disappearance of hysteresis), as Figll4l 
(b) shows, manages to be verified within the practically 
observable MC steps. In contrast, in the case of quasi 
2D systems expressed by ea. (|5n|l . it is quite difficult to 
verify such a full relaxation within practically possible 
MC steps because of a much larger number of degrees of 
freedom in a quasi 2D vortex state. 

Summarizing this subsection, for any |/3|-value with 
Hm lying below Hc2^ no true FOT occurs at Hc2, and 
the phase diagram Fig|Sl (a) is justified. For larger 
values when H„i may not be distinguished from Hc2^ it 
is practically difhcult to verify whether Fig|Hl(b), in place 
of FiglHl (a), is realized as a true phase diagram, and a 
hysteresis can appear in numerical experiments at H(,2 
even without a genuine FOT occurring there. 
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FIG. 14: (a) Data similar to FigHni (a) for the ID 

GL model, eq.lISnj, taken at MC steps 1 x 10'^. (b) Data 
at a — 4.74 showing a recovery of thermal equilibrium after 
many MC steps. 



IV. DISCUSSION 

As explained in Introduction, the present work was 
originally motivated as an extension of the problem of 
vortex phase diagram to the more general cases with 
spin (paramagnetic) depairing. Since the absence of MF 
second order transition at Hci in lower fields is well- 
established, it is unreasonable to expect the MF-FOT at 
i/c2 resulting from the Pauli paramagnetic depairing in 
higher fields to truly occur as a genuine FOT. The recent 
finding of an FOT-like nearly discontinuous crossover at 
i/c2 in the heavy fermion superconductor CeCoIns in 
H II c 0, provides us with a good occasion of a 
detailed comparison between the present theory and real 
data. Further, recent data in H _L c 0| showing a small 
hysteresis in heat capacity and suggesting a second or- 
der transition between the FFLO-like vortex solid and 
the ordinary vortex solid will stimulate detailed theoret- 
ical studies of vortex phase diagram in the region with 
paramagnetic depairing. 

As shown in Sec. II, in quasi 2D systems under fields 
perpendicular to the layers, a FOT at 7fc2 should nat- 
urally occur in T < T* at the mean-field level in clean 



enough superconductors with a moderately strong para- 
magnetic depairing, and the onset temperature Tfflo 
from which an FFLO-like vortex solid begins to appear 
lies much below T*. The results, that this modulated 
solid state is rarely seen compared with the apparent 
FOT at Hc2 and is easily lost by a small amount of 
impurities, are consistent with the observations |l3l | in 
a heavy fermion superconductor CeCoIns in H || c but 
seem to be a contrast to the conclusions on phase dia- 
gram in previous works Further, the range of T* 
(0.2 < T*/Tco < 0.36) shown in Fig.6 following from the 
value fioH°^^{0)/2TrTco = 0.8 and small (tTcq)" ^-values 
is comparable with th at { T* /Tro ~ 0.25) in CeCoIns 
under H || c [H [lljl6|. When the 7V(0)-value in 
CeCoIns in zero field [31| is used further, wc find the 
nearly discontinuous jump value of magnetization at Hc2 
when T = 50 (mK) to become 20 (G) which is compa- 
rable with the estimated value 30 (G) in ref.0|. Fur- 
ther, the width AH{T) of the magnetization jump at 
Hc2 will be estimated using the numerical data of Fie llll 
as follows. Using the relations ao{T) ~ 50(0.1Tco/T)2/3 
and P{T) cx T~^^^, suggested in Sec. Ill, and assuming 
|/3(T = 0.1Tco)| = 2.0, the data in Fig|TT]imply AH{T = 
O.lTco) ~ O.l(T) and AH{T = 0.025Tco) ~ 0.03 (T), 
which seem to be rather narrower than those of avail- 
able data 0,0|. It is not surprising because, even in 
a genuine FOT, an entropy (or a magnetization) jump is 
broadened due to an inhomogeneity in real materials. 

As shown and mentioned in Sec. Ill, the absence of a 
true FOT at Hc2 should be an appropriate interpreta- 
tion for observations 0,^^ at least in H || c where no 
measurable hysteresis was observed although the ordi- 
nary magnetic hysteresis related to the peak effect near 
Hc2 might appear. On the other hand, observa- 
tions of a t iny hysteresis were recently reported in s pec ific 
heat data [Ij, |l3] and also in magnetization data [I3 of 
CeCoIns in H _L c. Examining microscopic aspects in 
H _L c leading to a MF phase diagram is beyond the 
scope of the present paper and will not be considered 
here. However, the issue of the genuine phase diagram 
in Sec. Ill is more generic and may be applicable to the 
H _L c case in CeCoIus. Because the fluctuation effect 
in H _L c is weaker than that in H j | c at the same T, 
the observed hysteresis in refs.0,0,0| can be under- 
stood within the present theory arguing the absence of a 
genuine FOT at Hc2, if it has the same origin as a hys- 
teresis in numerical simulations at low enough T, shown 
in Figs ll3l and 1141 which arises from an incomplete re- 
laxation at long but finite time scales in a system with 
a strong MF-FOT. In relation to this, we point out that 
the onset of hysteresis in ref.|^lj| lies slightly above the 
temperature Tfflo at which i?FFLo(T) branches from 
the MF-FOT (i.e., Hc2{T)) Hue. It implies that the on- 
set of hysteresis has nothing to do with the appearance 
of the FFLO-like state and thus that there is no physical 
reason favoring a termination of some genuine SC (super- 
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conducting) FOT. Alternatively, the observed hysteresis 
might accompany a true magnetic FOT induced by the 
(|Ap) -nucreation at 7?c2 as a consequence of a coupling 
between the SC and a magnetic order parameters. In 
relation to this, we point out the observation of a large 
magnetic hysteresis [iTj in H _L c just below Hc2 to which 
the corresponding one was not seen in the vanishing of 
resistivity (i.e., a SC ordering). 

Finally, let us point out that the present theory easily 
explains why, in fields parallel to the layers, the transi- 
tion to an FFLO-like phase and the nearly discontinuous 
crossover at Hc2 implying the MF-FOT were measured 
not in organic materials with larger anisotropy but in 
a heavy fermion material with weaker anisotropy. At 
least at the MF level, the case with a field parallel to 
the layers in more anisotropic materials is subject to a 
stronger paramagnetic depairing and is the best candi- 
date for observing the paramagnetic effects such as the 
FFLO-hke state and the MF-FOT. The organic materials 
satisfy this requirement, and actually, the upwardly in- 
creasing Hc2{T) curve determined resistiviely ISfil u nder 
high fields parallel to the layers is an evidence [33 that 
the spin depairing is microscopically important without 
being disturbing by an impurity effect. However, both 
the MF-FOT and a transition to an FFLO-like state have 
not been seen in the organic materials. On the other 
hand, since such a upward ffc2(T)-curve is not visible in 
the heavy fermion material CeCoIns with a much weaker 
anisotropy, one might wonder why these crossover and 
transition arising from the spin depairing have occurred 
in this material. This puzzling facts are easily resolved 
by taking account of fluctuation effects examined in this 

Mer. Typically, in the organic and cuprate materials 
IssI , the fluctuation effect is much stronger compared 
with those of CeCoIns. Actually, a shorter coherence 
length tends to result in a larger Maki parameter and si- 
multaneously to enhance the fluctuation even in the par- 
allel field case 0,^9]. Consequently, as shown in Sec. Ill, 
the MF-FOT behavior is rounded and is transmuted into 
a broad continuous crossover by the SC fluctuation, re- 
flecting the absence of the true FOT at Hc2- Further, 
a remarkable field and temperature range of the vortex 
liquid region in which the resistance is finite may be cre- 
ated below Hc2{T)-cm:ve even in the parallel field case 
[3^ where the fluctuation effect is minimized. Since the 
FFLO-like state is limited to a narrow fleld range be- 
low Hc2, and the modulation parallel to the field does 
not lead to any ordering in the vortex liquid, the vor- 
tex liquid region should mask and erase the FFLO-like 
phase in a strongly fluctuating superconductor. For these 
reasons, cleaner superconducting materials with weaker 
fluctuation such as CeCoIn5 are the best candidates for 
examining the MF high field phase diagram in the case 
with Pauli paramagnetic effects. 
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APPENDIX A: DERIVATION OF K2 

In this appendix we present how to solve the eigenvalue 
problem of K2 or equivalently of D. Using the identity 
1/a = dpe~°'^ and after the energy integration, we 
get the differential operator of infinite order 



M2s) = J2 /°°dpe-(^^+V^)''Jo(2Jsin(^)p 



e>0 

X cos 



(2/p) (hppe— 



(54) 



where j7o is the zeroth order Bessel function. Expanding 
the exponential and averaging on the Fermi surface, we 
have 



(m!)2 + - 

+ off diagonal terms, (55) 

where 7r± are given below the eq. (17), and a circular 
Fermi surface was assumed. As explained above eq.(9), 
we have only to focus on the diagonal terms. Noticing the 

TV! 



(N- 



eigenvalue of 7r™7r™ in the N-th Landau level is 
and performing the m-summation, we finally obtain eq. 
(9). 



APPENDIX B: EXPRESSIONS OF h AND h 

In this appendix, we derive the expressions of I4 and /g 
in terms of a convenient expression for e^'^'^'^Uo,k{^±)- If 
we denote the position on a (2-dimensional) Fermi surface 
by a complex number vfC = I'a; 
pC* /th, we have 



iVy and define A 



-lA^. ^A7r+ ^A*Tr- 



(56) 



where we used the operator identity e^^^ = 
g- 2 [^'^Ig-^e^ valid when is a constant. From 

this expression it is sufficient to know the form of 



e%/2'^'^+UQ^(rj^). Using the above identity to obtain 



(57) 
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and noticing e"^-g{x) = g{x + a) for any non-singular e'''^ "*Uo ^(r_L) = e^^^l^l'"^* )e-5 W'-ff+fc'-H-A*) -iky ^ 
function we finally have j^gg-j 

(58) 



With the help of the above identities, the following results are easily derived. 
^-L i(piv-nT+P2v-n2+P3v-n.^+P4v-n4) * 



S 



H 



^wS,fci(l-±l)wo,/c2(l-J^2)M5,fc3('"J-3)"0,fc4(l"±4) 



1 

1 

71 



Z]t=l l^»l^ + 5(^l3+^24) + (^t+\D(-^2 + A4)+2rff (fcl3A*3-/C24A24) 



(60) 



Sh 



W0,fci(l"±l)w0,fc2(r-L2)uS,fe3(r±3)M0,/c4(l"±4)uS,fc5(l"-L5)"0,fc6(l"±6) 



%/3 



<5/ci+fc3 + fc5:fe2+fc4+fc6S 



xe 



r5^i = ll-*>i|^ + i(i;^dd-'^i^+i:ovcn-'*i )^ A(i:^dd "^i + i :cvon '''' ) ^ ~ I ( (i , jjJodd (i ,j)^evei, "''ij ) "I f"((i,j)^odd '^ij •''ij ~ (i , j)^ovon '='3 J ) 



1 -—V 

= ;y|4i+fc3+fc5,fc2+'c4+fc6e ' ^<''^' ''^({Aj}) 
where (», j) = {(1, 3), (3, 5), (5, 1), (2, 4), (4, 6), (6, 2)}. 



(61) 
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